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The densest amorphous packing of rigid particles is known as random close packing. It has long 
been appreciated that higher densities are achieved by using collections of particles with a variety of 
sizes. The variety of sizes is often quantified by the polydispersity of the particle size distribution: 
the standard deviation of the radius divided by the mean radius. Several prior studies quantified the 
increase of the packing density as a function of polydispersity. Of course, a particle size distribution 
is also characterized by its skewness, kurtosis, and higher moments, but the influence of these param- 
eters has not been carefully quantified before. In this work, we numerically generate many packings 
with different particle radii distributions, varying polydispersity and skewness independently of one 
another. We find two significant results. First, the skewness can have a significant effect on the 
packing density and in some cases can have a larger effect than polydispersity. Second, the packing 
fraction is relatively insensitive to the value of the kurtosis. We present a simple empirical formula 
for the value of the random close packing density as a function of polydispersity and skewness. 



I. INTRODUCTION 

Understanding various aspects of random close pack- 
ing (rep) has great scientific and industrial importance 
[1] as it has been linked to a wide range of problems such 
as the structure of living cells [2], liquids [31 0], gran- 
ular media [SHE]) emulsions [5], glasses [TU], amorphous 
solids [TT], jamming [TH [T3], and the processing of ce- 
ramic materials |14j . Random close packing is typically 
defined as a collection of particles packed into the dens- 
est possible configuration, although more rigorous defini- 
tions are available [15] . Experiments have found that the 
densest random packing of monodisperse spheres typi- 
cally occurs close to <po, rcp ~ 0.64 [JJ, where the density 
4> (or packing fraction) is defined as the ratio of the total 
volume occupied by spheres to the volume of the con- 
tainer. 

Formally, a packing consists of particles with a dis- 
tribution in radii P(R). The polydispersity is defined 
as 6 = y/(AR 2 )/(R). Here, AR = R - (R) and the mo- 
ments of R (and AR) are defined as (R n ) = J R n P(R)dR 
(and (AR n ) = / R n P(R)dR). It has long been appreci- 
ated that packings of spheres can have other rep densities 
when 5 > [T6Tf2"T] . Prior experiments 22 26 and simu- 
lations [2TH5T] have nicely shown that as the polydisper- 
sity increases, the particles pack to higher volume frac- 
tions because the smaller particles pack more efficiently 
by either layering against larger particles or by fitting 
into the voids created between neighboring large parti- 
cles [231 I32"H34"] . In practice, depending on the degree of 
the polydispersity, the packing fraction can increase from 
0.64 for monodisperse packings to nearly ~ 0.75 for pack- 
ings with 0.65 polydispersity [53]. For the extreme case 
of two different particle sizes with a size ratio approach- 
ing infinity, the voids between the large particles can be 
packed randomly with small particles and so 4> can be as 
large as o ,rc P + (1 - 0o,rc P )0o,rc P ~ 0.88 [231[35]. 

While it is intuitive that the polydispersity can af- 
fect 4> rcp , it is also reasonable that the shape, not just 
the spread, of the distribution P(R) may also influence 



<p r cp L l, 23J . For instance, an infinite number of distribu- 
tions can have the same value of 8 but yet differ in their 
form. One can characterize the shape using the skewness 
S = (AR 3 )/(AR 2 ) 3 / 2 , kurtosis K = (AR 4 ) / (AR 2 ) 2 , and 
higher moments. There have been prior studies that have 
investigated the influence of distribution shape on the 
density of tightly packed particles QH H3 OS I35PD]. 
Similar to the studies on polydispersity, they find that 
the shape of the particle distribution can have a pro- 
found influence on the packing density. However, these 
prior studies either did not independently vary S and S 
but rather conflated the influences of both, or else used 
other metrics besides S and S to quantify P(R). Of the 
prior studies, Tickell et al. [UJ is the only one to report 
on the effects of skewness and kurtosis for experiments 
carried out with sand, finding that over a narrow range 
in skewness the packing density can increase by 0.04 with 
no dependence on kurtosis. However, they did not con- 
trol for polydispersity, leaving it unclear the relative im- 
portance of polydispersity and skewness. The key unan- 
swered question by the prior work is how the skewness 
of a distribution influences 4> rC p, and large this effect is 
relative to the effects of polydispersity. 

To study the influence of skewness on <p rcp , we generate 
random close packings numerically, and control for poly- 
dispersity and skewness independently to address how 
dominant each term is in determining 4> rcp . Our method 
for generating these packings was previously developed 
in Ref. [35]. Briefly, infinitesimal particles are placed 
randomly in a periodic container, gradually expanded, 
and moved at each step to prevent particles from over- 
lapping. At the beginning of the simulation, particles 
are assigned radii with a specific distribution and as the 
particles expand they do so by a multiplicative factor 
such that the shape of the radii distribution is fixed. The 
value of 4> rcp is known to be sensitive to protocol [J3J [33] , 
and it is not known if this algorithm or any other algo- 
rithm produces rigorously defined random close packed 
states [HI 221 [45l 06]. Our goal is not to determine 
the precise value of (prep for a given P(R), but rather 
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FIG. 1. (Color online), (a) Examples of three different par- 
ticle radii distributions with the same polydispersity of 0.25 
and nearly the positive same skewness. The binary distrub- 
tion has S = 0.78, the linear distribution has S = 0.56, and 
the lognormal distribution has S = 0.78. (b) Examples of 
two different particle radii distributions with the same poly- 
dispersity of 0.25 and negative skewness S = —0.5. 



to empirically understand the trend in <f> rcp with poly- 
dispersity and skewness. Nonetheless, we note that our 
algorithm when applied to a monodisperse packing gives 
4>o, rep ~ 0.64, close to the experimentally found value 
and in agreement with prior simulation work. 

To efficiently determine 4> rC p for a chosen particle size 
distribution, we exploit the known finite size dependence 
4>rc P (h) = (pr^p — C/h, where h is the system size, is 
the random close packing fraction in the limit h — > oo, 
and C is a fitting constant [17] • By generating many 
packings with different box sizes h, we fit 4> rC p(h) to de- 
termine for each distribution. We generate packings 
with box sizes of ~10, 14, 18, and 23 mean particle di- 
ameters in length to determine (f)^. For the rest of the 
paper, 4> rcp will be used to indicate (f>^ p . 

To control for both 5 and 5 independently, we study 
packings using four different distributions: binary, linear, 
gaussian, and lognormal. The binary and linear distribu- 
tions are determined by two control parameters, allow- 
ing for us to control S and 5 independently, while the 
gaussian and lognormal distributions are determined by 
only one parameter, and therefore 6 and 5 can not be 
controlled independently. By generating many packings 
with different 6 and 5 using these four distributions, we 
can compare the results to see how sensitive 4> rcp is to 
polydispersity and skewness, but we can also compare 
different distributions with the same 8 and 5 to see how 
sensitive <j) rcp is to other subtle differences in the distri- 
bution shape. For all distributions, we impose (R) = 1. 

More specifically, the binary distribution consists of 
particles with two distinct radii. The shape of the distri- 
bution is determined by the size ratio and number ratio of 
these two particle types. The linear distribution is a con- 
tinuous distribution of the form P(R) = AR + B, where 
the distribution in particle size exists over a finite range 
a < R < b. Our choice of (R) = 1 and the requirement of 
normalization (J P(R)dR — 1) imposes two constraints 
on the parameters (a, b, A, B). For the two remaining de- 
grees of freedom, we define 77 = b/a and p = P(a)/P(b). 



We compute 5 and 5 for a grid of r\ and p values, and 
then interpolate to find the parameters for P(R) for the 
desired 5 and 5 values, allowing us to vary them system- 
atically. The third distribution is a Gaussian of the form 
P(R) = Aq exp(— [R— 1) 2 /2ct 2 ), where a is the standard 
deviation and Aq = 1/ \<j\/2k). For larger a, some of the 
particle radii could be negative, which is unphysical, or 
very close to zero, which may prevent generating packings 
within a reasonable timeframe. To avoid these issues, we 
truncate the Gaussian distribution such that the small- 
est particle radius is no smaller than 0.1. The Gaussian 
distribution has a fixed skewness 5 = except for the 
trunctated Gaussians, which have a slight positive skew- 
ness. The last distribution we consider is the lognormal 
distribution P(R) = A L exp(-0.5(ln R/a + 0.5a) 2 )/R, 
where Al = l/(a^/2ir). Similar to the Gaussian distri- 
bution, the skewness of the lognormal distribution is not 
adjustable, but is always positive and becomes larger as 
a becomes larger. 

In Fig. [lja), we compare three different distributions 
with polydispersity 5 = 0.25 and nearly the same pos- 
itive skewness 5 0.75. We see that the distributions 
are quite different, in particular in their tails. For exam- 
ple, the linear distribution has many more small particles 
than the other two distributions. The lognormal distribu- 
tion has tails that include both smaller and larger parti- 
cles than the other two distributions. It's not necessarily 
obvious how the values of <j) rC p will be ranked for these 
cases. In Fig. Wb), we show two different distributions 
with polydispersity 6 — 0.25 and skewness 5 rs —0.5. As 
these distributions have negative skewness, both distri- 
butions have more larger particles than smaller particles. 
Once again, it's not necessarily clear how <p r cp should 
differ between the two packings. 

After generating nearly 10,000 packings with different 
particle radii distributions, we plot (f> rcp as a function of 
skewness 5 for all our data in Fig. [2] with the different 
groups of data (different colors) corresponding to differ- 
ent polydispersity values 5. Each data point in the figure 
has a one to one correspondence to the distribution type, 
5, and S. The symbol or line type of the data indicates 
the P(R) distribution type. Remarkably, the figure shows 
that regardless of the type of particle radii distribution, 
(f>rcp is nearly the same for the same pairing of polydis- 
persity and skewness. It also shows that tp rcp increases 
with both increasing S and S. Strikingly, the skewness 
can have an equally important effect as the polydisper- 
sity. For example, for S — 0.40 and 5 = 0, <fi r cp is shifted 
upward by ~ 0.02. Fixing that value of 5, changing 5 
to ±1 shifts tfircp by w ±0.02. For highly polydisperse 
samples, one cannot accurately know 4> rcp without also 
knowing the skewness of the radius distribution. For the 
binary samples (solid lines in Fig. |2| 5 can be even larger 
in magnitude and have an even larger influence on <j) rcp 
than 8 has. 

The increase in 4> rcp with skewness is not uniform. For 
negative skewness (more big particles), the polydisper- 
sity 5 does not seem to influence <j) rC p as much as when 
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FIG. 2. (Color online). This figure shows how <f> rcp depends 
on particle size distribution, polydispersity S, and skewness 
S. The solid lines represent <f) rcp for binary packings and 
the symbols represent (f> rcp for packings with either linear, 
Gaussian, or lognormal particle distributions as indicated by 
the legend. The colors represent different polydispersities of 
either 0.1, 0.15, 0.25, or 0.4. The dashed lines are a fit to the 

= 0.634, 



data using 



+ ci<5 + C2SS , where 



d = 0.0658, and c 2 = 0.857. 



the skewness is positive. This is not too surprising since 
the volume of each particle grows with R 3 . When the 
total number of bigger particles is greater than the to- 
tal number of smaller particles (negative skewness), the 
volume occupied by all the large particles is significantly 
greater than the volume occupied by all the small parti- 
cles. In effect, the big particles pack like a low polydis- 
persity sample and occupy the majority of the container, 
while the small spheres occupy an insignificant portion, 



and 



J rcp 



approaches 



0.64 for a monodisperse 



sample. For positive skewness (more smaller particles), 
<p rc p has a fairly strong dependence on 5 and S, where 
4> rC p increases with increasing number of small particles. 
The reason for this increase in <p rcp is likely due to the 
small particles fitting into the spaces between larger par- 
ticles. As discussed in prior work [231 133 03] > the lo- 
cal porosity is smaller around two neighboring particles 
of different sizes than around two neighboring particles 
of the same size. This effect is greater for larger dif- 
ferences in the size of two neighbors. As skewness and 
polydispersity increase, both the number of small par- 
ticles present and the average size discrepancy between 
neighboring particles increase, resulting in a larger 4> rcp . 

To provide a qualitative sense of the behavior, Fig. [3j 
shows a 2D slice through four different 3D packings, (a) 
and (b) are lognormal packings, where (a) is a packing at 
low polydispersity and skewness and (b) is a denser pack- 
ing at a higher polydispersity and skewness. Packings (c) 
and (d) are two binary packings with polydispersity 0.4, 
where (c) has a large negative skewness and (d) is denser 
and has a large positive skewness. As discussed above, 
at large 5 and S, small particles can either layer around 
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FIG. 3. (Color online). Each image represents a 2D slice 
through a 3D packing, and the dashed box is the boundary of 
the periodic packing. The volume fraction for each packing 
shown is close to the extrapolated 4>rc P - Also, the area fraction 
of each 2D slice is the same as the volume fraction of the 3D 
packing they represent. 



larger particles and/or fit in the voids between bigger 
particles. In Fig. [3^b) and (d) we see evidence of small 
particles sitting in the void areas between big particles. 
In Fig. j3j[a) and (c), where the skewness is lower, we see 
less evidence of this. These observations are consistent 
with the results of Fig. [2j 

Since 4> rcp is nearly determined by the two parameters 
6 and S, we fit all the data to a simple equation (f> rc p = 
4>* cp + c%S + C2SS 2 , where (f>* is the packing fraction for 
a monodisperse packing of spheres (5 = and S = Q) 
and Ci and C2 are empirical constants. These fit lines are 
shown as dashed lines in Fig. [2j and agree reasonably well 
with the data. We also tried fits with higher order terms 
in S and S, but we found first order in S and second order 
in 6 reasonably fit all the data well. Using our simple fit 



model, we find </>* cp = 0.634, c x = 0.0658, and c 2 = 0.857. 
The fit gave a value for (f>* cp close to the experimentally 
accepted value of 0.637 [Tl kE8"]. 

There arc slight differences in the </> rcp values for differ- 
ent distribution types for the same 5 and S values, seen 
in Fig. [2j Thus far we have focused on 5 and S to char- 
acterize our distributions, and of course the distributions 
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differ in their higher moments. The next quantity to con- 
sider is the kurtosis K defined above, and it might be a 
potentiai additional parameter to explain the variations 
in Fig. [2j However, comparing packings with the same 
8 and S we find that the slight differences in (j) rcp have 
no clear dependence on K, suggesting that the kurtosis 
is not an important contributor to the <f> rcp . This agrees 
with the 1933 qualitative observations of Tickell et al. 

0H- 

We note in passing that a rapid prediction of <j> rcp can 
be found by using the numerical algorithm of Farr and 
Groot [35], based on any P(R) as input. All of our nu- 
merical data from explicitly generating 3D packings agree 
quite well with the results of their algorithm. 

Our data have two significant conclusions. First, the 
skewness S has a significant influence on 4> rcp for dis- 
tributions with a large polydispersity <5. Second, know- 
ing S and S can determine <j) rcp to within approximately 
±0.002 without taking into account any other details of 
the shape of P(R). 

This collapse of 4> rcp values for a given 6 and S but 
different distribution shapes is intriguing, as presumably 



the structures within the packings are different for dif- 
ferent P(R). For that matter, one can have the same 
4> rcp value for different 6 and S, see for example Fig. ^ a) 
and (c), and clearly these will have different microstruc- 
tures. This might be useful for studying aspects of the 
jamming transition of spherical particles. Many prior 
results show that various properties of these systems 
depend on the distance to the jamming point |49H53j . 
where the jamming point is thought to be the same as 
4>rc P [HI H3 [46] [54] . One can imagine conducting ex- 
periments or simulations to compare the properties of 
packings near the jamming transition with different mi- 
crostructures, but the same jamming point. These could 
be equally useful for studying the colloidal glass transi- 
tion, which may be influenced by 4> rcp [5TJ I55H59] . Such 
experiments may provide further insight into the univer- 
sal nature of the jamming transition and glass transi- 
tion, but may also highlight subtle dependencies on the 
microstructure . 
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